# -*- coding: utf-8 -*-
"""
Created on Mon Jan 11 20:17:59 2021

@author: Lenovo
"""

import numpy as np
from numpy import random
import matplotlib.pyplot as plt
global Gm,F0,Np,CR,D
import math


#加载数据集到dataMat,跳过第一行标题，
dataMat=np.loadtxt(open('D:\\冠心病稳定期VD.csv'),delimiter=',',skiprows = 1)
row = dataMat.shape[0] #数据集行数，数据条数
col = dataMat.shape[1] #数据集列数，症状个数
#############参数
F0 = 0.5 #变异概率
CR = 0.2 #交叉概率
G = 1    #当前进化代数
D = 2*col + 1   #个体维度、等于2倍症状个数加1
Np = 8*D  #种群规模
#下面参数可调整

Gm=200#最大进化代数
E1=0.5#主症的阈值
E2=0.1#次症的阈值

V=np.zeros((Np,D))
U=np.zeros((Np,D))
XG_next=np.zeros((Np,D))
value=np.ones(Np)*-float('inf')#初始化都为负无穷
cv=np.zeros(Np)#存储罚分值，该值越小越好
cv[:]=1000#初始化罚分值都为1000
PF=np.zeros(Gm)

##############产生初始种群##
max0=1
min0=0
X0 =random.random(size=(Np,D))
#print(X0)  #随机初始化v0-v16，他们的取值范围都为[0-1],也可以写到一个循环里。
#需要注意的是每连续的两个是一个症状，例如v0和v1是一个症状、v2和v3是一个症状，最后v26是类别概率值
for i in range(27):  #x0--x27
    X0[:,i] = X0[:,i]*(max0-min0) + min0
XG=X0
#print(XG)

###################计算函数值，返回每个个体的似然度值
#v[j]中，范围是0-26，每个症状由2个参数，每连续的2个值是属于同一个症状的，除最后一个是类别参数
# j是偶数时p(x=1|z=1)概率；j是奇数时为p(x=1|z=0)的概率，
#例如v[0]是p(症状x0=1|z=1)概率，v[1]是p(症状x0=1|z=0)概率，
#  v[2]是p(症状x1=1|z=1)概率，v[3]是p(症状x1=1|z=0)概率，
#  v[4]是p(症状x2=1|z=1)概率，v[5]是p(症状x2=1|z=0)概率，
def fcomput(v,index):
    allsumlikelood = 0 #个体的BIC分，等于所有条数据的似然度之和再减去罚分
    sumlikelood = [] #个体在每条数据上的似然度
    #penalize是计算BIC分需要的后半部分，参数个数D除以2，然后乘以log(row)
    penalize = (D/2)*(math.log(row)) #对模型复杂程度的罚分
    for i in range(row): #对数据集中的每条数据
        #先算p(z=1)情况下概率
        likeloodforZis1 = 1
        likeloodforZis1 = likeloodforZis1 * v[16]
        for j in range(col):
            if dataMat[i][j] ==1:#p(x=1|z=1)时
                likeloodforZis1 = likeloodforZis1 * v[2*j]
            else: #p(x=0|z=1)时
                likeloodforZis1 = likeloodforZis1 * (1-v[2*j])
        #print(likeloodforZis1)
        #先算p(z=0)情况下概率 
        likeloodforZis0 = 1
        likeloodforZis0 = likeloodforZis0 * (1-v[16])
        for j in range(col):
            if dataMat[i][j] ==1:#p(x=1|z=0)时
                likeloodforZis0 = likeloodforZis0 * v[2*j+1]
            else: #p(x=0|z=0)时
                likeloodforZis0 = likeloodforZis0 * (1-v[2*j+1])# #p(x=0|z=0)
        sumlikelood.append(likeloodforZis0+likeloodforZis1)
    for i in range(row):
        allsumlikelood += math.log(sumlikelood[i])
    print('第%d个个体的似然度是' %index,end='')
    print(allsumlikelood-penalize) #似然度之和再减去罚分，即是BIC分
    return (allsumlikelood-penalize)

########计算约束惩罚度
def g(v):
    g1=((v[0]-v[1])/(1-v[1]))-E1#5个主症，满足每个主症的公式大于阈值
    g10=((v[18]-v[19])/(1-v[19]))-E1
    g11=((v[20]-v[21])/(1-v[21]))-E1
    g12=((v[22]-v[23])/(1-v[23]))-E1
    g13=((v[24]-v[25])/(1-v[25]))-E1
    
    g2=E2-((v[2]-v[3])/(1-v[3]))#8个次症，满足每个主症的公式小于阈值
    g3=E2-((v[4]-v[5])/(1-v[5]))
    g4=E2-((v[6]-v[7])/(1-v[7]))
    g5=E2-((v[8]-v[9])/(1-v[9]))
    g6=E2-((v[10]-v[11])/(1-v[11]))
    g7=E2-((v[12]-v[13])/(1-v[13]))
    g8=E2-((v[14]-v[15])/(1-v[15]))
    g9=E2-((v[16]-v[17])/(1-v[17]))
    
    g1=min(g1,0)#5个主症
    g10=min(g10,0)
    g11=min(g11,0)
    g12=min(g12,0)
    g13=min(g13,0)
    
    g2=min(g2,0)#8个次症
    g3=min(g3,0)
    g4=min(g4,0)
    g5=min(g5,0)
    g6=min(g6,0)
    g7=min(g7,0)
    g8=min(g8,0)
    g9=min(g9,0)
    g=min(g1,g2,g3,g4,g5,g6,g7,g8,g9,g10,g11,g12,g13)
    if g==0:
        cv=0
    else:
        cv=pow(g1,2)+pow(g2,2)+pow(g3,2)+pow(g4,2)+pow(g5,2)+pow(g6,2)+pow(g7,2)+pow(g8,2)\
            +pow(g9,2)+pow(g10,2)+pow(g11,2)+pow(g12,2)+pow(g13,2)
    return cv

#进化开始
while G<=Gm:
    ####变异操作部分#######################
    #print(G)
    for i in range(Np): 
        li= [] #list(range(Np))
        for temp0 in range(Np):
            li.append(temp0)
        random.shuffle(li)#打乱个体的序号
        #print(li)
        dx=li
        j=dx[0]
        k=dx[1]
        p=dx[2]
        #要保证与i不同
        if j==i:
            j=dx[3]
        elif k==i:
            k=dx[3]
        elif p==i:
            p=dx[3]
        #具体变异操作
        suanzi= math.exp(1-float(Gm)/(Gm+1-G))
        F=F0*(2**suanzi)
        mutant=XG[p]+F*(XG[j]-XG[k])
        for j in range(D):
            if mutant[j]>min0 and mutant[j]<max0:
                V[i,j]=mutant[j]
            else:
                V[i,j]=(max0-min0)*random.rand()+min0

    #####################下面开始交叉操作#####################   
    for i in range(Np):
        randx= [] #list
        for temp in range(0,D-1,2):#D=27维,，
            randx.append(temp)
        random.shuffle(randx)#随机打乱randx即维度D的序号
        #D=27维,j范围是0-26，隔一个取一个，这很重要，
        #因为1个症状有两个参数，即每连续的2个值是属于同一个症状的。除最后一个是类别
        for j in range(0,D-1,2):
            if random.rand()>CR and randx[0]!=j:
                U[i,j]=XG[i,j]   #同一个症状的两个参数要同时选过去
                U[i,j+1]=XG[i,j+1]
            else:     #同一个症状的两个参数要同时选过去
                U[i,j]=V[i,j]
                U[i,j+1]=V[i,j+1]
        #下面是对最后一个j=26时，即类别p(z)的值进行交叉
        if random.rand()> CR:
            U[i,26]=XG[i,26]
        else:    
            U[i,26]=V[i,26]
    ###################下面开始选择操作#############
    for i in range(Np):
        f1=fcomput(XG[i],i)
        f2=fcomput(U[i],i)
        g1=g(XG[i])
        g2=g(U[i])
        if g1==0 and g2==0:
            if f1>f2:
                XG_next[i]=XG[i]
            else:
                XG_next[i]=U[i]
        elif g1==0 or g2==0:
            if g1==0:
                XG_next[i]=XG[i]
            else:
                XG_next[i]=U[i]
        elif g1!=0 and g2!=0:
            if g1<g2:
                XG_next[i]=XG[i]
            else:
                XG_next[i]=U[i]   
    ############下面更新代数，往前进化1代
    for i in range (Np):
        value[i]=fcomput(XG[i],i)
        best_value=min(value)
    PF[G-1]=best_value
    XG=XG_next  
    G=G+1
####进化结束
with open('DE.txt','w') as f:
    for i in PF:
        print(i, file=f)
plt.plot(PF)
plt.show()
print('进化结束')
for i in range(Np):
    cv[i]=g(XG[i])
    if abs(cv[i]-0)<0.005:
        value[i]=fcomput(XG[i],i)
        #print('cv[i]为0时',end='')
        #print(value[i])
best_value=max(value)
temp=value.tolist()
pos_max=temp.index(max(temp))
print ('求得f似然度最大(BIC)为: %f'%best_value)
print('f最大时解的违约罚分(最小)为：%f' %min(cv))
best_vector=XG[pos_max]
print ('f最大值时对应的解为：',end='')
for i in range(len(best_vector)):
    print('%0.2f, ' %best_vector[i],end=' ')#小数点后保留2位
print('')


########################################下面计算阈值和分值
###########下面计算阈值
thefirst = math.log((1-best_vector[D-1])/best_vector[D-1]) #计算阈值前半部分
sumv = 0 #计算阈值公式中的后半部分
for i in range(0,D-1,2):
    temp = (1-best_vector[i])/(1-best_vector[i+1])
    sumv = sumv + math.log(temp)
thesorld = thefirst - sumv 
scalerate = thesorld/10 #缩放率等于原始阈值除以10，这样把阈值统一为10        
###########下面计算每个症状的分值
name=['胸闷','肢体困重','体胖','口黏','面色浑浊','大便黏滞','脘腹痞满','嗜睡','纳差','舌胖边有齿痕','苔腻','苔白滑','脉濡或滑']
score = 0
for i in range(0,D-1,2):
    s = (best_vector[i]/(1-best_vector[i]))/(best_vector[i+1]/(1-best_vector[i+1]))
    score = (math.log(s))/scalerate #原始分值除以缩放率，这样把阈值统一为10
    j = (int)(i/2) #症状有8个，变换为症状的下标
    print('%s:' %(name[j]), end='')
    print('%0.2f' %score) #小数点后保留2位